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In Phys. Rev. Lett. 114 (2015) 6, 061601, we reported on a new universality class for longitudinally 
expanding systems, encompassing strongly correlated non-Abelian plasmas and iV-component self¬ 
interacting scalar field theories. Using classical-statistical methods, we showed that these systems 
share the same self-similar scaling properties for a wide range of momenta in a limit where particles 
are weakly coupled but their occupancy is high. Here we significantly expand on our previous work 
and delineate two further self-similar regimes. One of these occurs in the deep infrared (IR) regime of 
very high occupancies, where the nonequilibrium dynamics leads to the formation of a Bose-Einstein 
condensate. The universal IR scaling exponents and the spectral index characterizing the isotropic 
IR distributions are described by an effective theory derived from a systematic large-expansion at 
next-to-leading order. Remarkably, this effective theory can be cast as a vertex-resummed kinetic 
theory. The other novel self-similar regime occurs close to the hard physical scale of the theory, 
and sets in only at later times. We argue that the important role of the infrared dynamics ensures 
that key features of our results for scalar and gauge theories cannot be reproduced consistently in 
conventional kinetic theory frameworks. 
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I. MOTIVATION AND OVERVIEW 

Universality far from equilibrium. In a recent letter we 
demonstrated universal behavior, in a wide range of mo¬ 
menta, for self-interacting scalars and non-Abelian plas¬ 
mas when they are longitudinally expanding [T]. This is 
remarkable because the scattering processes controlling 
the dynamics of scalar and gauge field theories are in 
principle very different from one another. The universal 
behavior is observed for systems with sufficiently high oc¬ 
cupancy / ^ 1 in the limit of weak coupling A ^ 1. Here 
A denotes the quartic self-interaction for scalar fields; for 
non-Abelian plasmas, it is replaced by the gauge coupling 
^ ag. We showed that the far-from-equilibrium dynam¬ 
ics of the different systems in the high occupancy limit 
is governed by a common nonthermal fixed point with 
the same dynamical scaling exponents that control the 
time evolution of characteristic longitudinal and trans¬ 
verse momentum scales. Even detailed properties of the 
different systems, such as the single particle distributions 
in the given momentum range, turn out to have identi¬ 
cal profiles. This opens up the possibility of studying 
important aspects of the dynamics of such complex sys¬ 
tems using the simplest representative field theory of the 
corresponding universality class in a well controlled weak- 
coupling limit. 

Classieal-statistieal field theory. A major motivation 
for these studies are experiments colliding heavy nuclei 
in the ultrarelativistic high energy limit. In the early 
stages of such a collision, highly occupied gluon fields 
with typical momentum Q are expected to be formed 0- 
HS]. Since the typical occupancies ^ l/as{Q) are large, 
the gauge fields can be strongly correlated even for small 


gauge coupling as{Q). In such highly occupied systems 
dynamical quantum corrections are suppressed at weak 
coupling and the nonequilibrium quantum dynamics in 
this transient regime can be accurately mapped onto a 
classical-statistical problenj^ The latter can be solved 
using real time lattice simulation techniques HiHaS]. 

The classical-statistical field description is limited to 
the early stages of the evolution since the typical occu¬ 
pancies decrease at later times. It is well known that 
classical dynamics cannot describe the late time evolu¬ 
tion towards thermal equilibrium due to the Rayleigh- 
Jeans divergence nor can it isotropize in the presence of 
rapid longitudinal expansion. The latter property has 
been confirmed in classical-statistical simulations of the 
Yang-Mills theory in 1211122 ] and as well as in classical 
kinetic theory simulation^ of gauge [38] and scalar the¬ 
ories [39] . 

Effeetive kinetie theory. Kinetic theory can in principle 
describe the late time behavior of systems if the typical 
occupancies are not too large (/ ^ 1/A) and if their typi¬ 
cal momenta are much larger than the in-medium screen¬ 
ing scale mo, such that ^ ^ f f{p)/P^ where 

denotes the phase space integration. While these con¬ 
ditions for kinetic theory are not fulfilled for / ^ 1/A at 


^ In this regime, quantum corrections are suppressed by one order 
of the coupling constant even in the anisotropically expanding 
theories. 

^ Stated differently, isotropization/thermalization from classical 
simulations introduces spurious UV divergencies. In that event, 
the results no longer describe the corresponding quantum field 
theory dynamics IM1I571. 
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early times, both kinetic theory and classical-statistical 
field theory can in principle give equally valid descrip¬ 
tions for a wide range of typical perturbative occupancies 
1</<1/A goHil]. 

A number of weak coupling thermalization scenarios 
of gauge theories, based on kinetic frameworks, have 
been developed over the years. These differ primarily in 
how infrared momentum modes below the scale ttid are 
treated jdSHiO] . leading to very different paths whereby 
thermalization is achieved. Since classical-statistical field 
theory accurately captures the physics of the nonpertur- 
bative infrared (IR) regime, it was successfully employed 
in ED ED to single out the kinetic description of [43] 
- henceforth referred to by the acronym BMSS. As we 
shall discuss shortly, this result is surprising and requires 
a deeper understanding of the underlying dynamics. 

We should nevertheless first observe that our results 
suggest that a comprehensive strategy for weak coupling 
thermalization is to start with classical-statistical field 
theory at early times and to continue later with an ap¬ 
propriate kinetic theory, passing from one regime to the 
other in the time window where both descriptions are 
valid. This program has led to significant progress in our 
understanding of the thermalization process of expand¬ 
ing non-Abelian plasmas as summarized in [SHllSollSI]. 
The effective kinetic theory description of [38] accu¬ 
rately reproduces the classical-statistical results in the 
expected regime for weak coupling and sufficiently large 
momenta. This concerns even detailed properties such as 
the scaling and shape of the single particle distribution in 
that range. If applied to stronger couplings, the kinetic 
descriptions show significant deviations from classical- 
statistical behavior because of the additional quantum 
corrections [38] [39] that now become of the same order as 
the classical contributions. Before extending the results 
to stronger couplings, it is crucial to fully understand the 
well controlled weak coupling regime. 

Puzzles. The result from classical-statistical simula¬ 
tions, demonstrating that BMSS provides the correct ki¬ 
netic description, was unanticipated. This is because a 
consistent kinetic description of expanding non-Abelian 
plasmas would suggest that plasma instabilities [52H56] 
should qualitatively modify the evolution even at late 
times nil Hi im El ED . However classical-statistical 
simulations of ED El suggest that plasma instabilities 
play no significant role beyond the very early stages. 
While a microscopic understanding of the absence of late 
time instabilities remains an open question, this obser¬ 
vation has been incorporated in the kinetic description 
of [38]. 

A further interesting puzzle for longitudinally expand¬ 
ing systems concerns the ratio of the longitudinal pres¬ 
sure to the transverse pressure, whose value is a 

measure of the degree of isotropization in the system. In 
kinetic theory, for non-Abelian gauge theories, small an¬ 
gle scattering dominates for proper times r < Q~^as 
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FIG. 1. Comparison of the evolution of the longitudinal to 
transverse pressure ratio in longitudinally expanding scalar 
and gauge theory. While parametric estimates suggest a scal¬ 
ing as Pl/Pt ~ in both cases, the observed decay 

of the longitudinal pressure is significantly slower due to large 
infrared contributions to the longitudinal pressure. This puz¬ 
zle is elaborated on further in the text. 


Parametric estimates a la BMSS would then give 


kinetic 

L theory 


(Qr)-2/3. 


( 1 ) 


However despite the fact that kinetic and classical- 
statistical descriptions seem to agree well for momenta 
above the screening scale Eaisn], the field theoretic sim¬ 
ulations of the non-Abelian plasma [H] [22] show signifi¬ 
cant deviations from the kinetic theory scaling of Eq. 0- 
This is demonstrated in Fig.[^ where the behavior of the 
gauge theory [22] is shown for isotropic initial conditions. 

In contrast to the gauge theory case, the matrix ele¬ 
ment for self-interacting scalars shows no preference for 
either small or large angle scattering. Despite this differ¬ 
ence, kinetic theory simulations of the expanding scalar 
theory [39] obtain exactly the result shown in Eq. 0 . 
This agreement with the gauge theory is surprising and 
not evident from analytical considerations so far. Be 
that as it may, the classical statistical simulations for 
the scalar theory shown in Fig. 0 deviate significantly 
from this kinetic theory result for large initial occupancy. 

Since the pressure is a bulk quantity, which sums over 
all momentum modes, an immediate conclusion from the 
observed discrepancy is that the infrared regime leads to 
significant contributions to Pl/Pt- The scalar example 
is useful because the influence of infrared dynamics on 
the evolution of bulk quantities can be studied without 
the gauge ambiguities of Yang-Mills theories. 

Infrared dynamies and Bose eondensation. There are 
a number of studies suggesting that there are important 
corrections from infrared modes or Bose condensates at 
transient stages of the thermalization process for non- 
Abelian plasmas as well [dS] [58ll64] . However a nonper- 
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turbative analysis explaining the classical-statistical sim¬ 
ulations results is challenging for the gauge theory. As 
noted, the situation is very different for scalar field theo¬ 
ries. In this case, a detailed analytical understanding can 
be achieved directly based on the underlying quantum 
field theory. The observed universality between longitu¬ 
dinally expanding scalars and non-Abelian plasmas [1] in 
a wide range of momenta could possibly be put to use 
to obtain further valuable insights into the gauge theory 
dynamics. 

One may also try and simulate infrared effects within 
kinetic frameworks for scalar field theories. One pop¬ 
ular kinetic description for scalar field theory, recently 
implemented for a longitudinally expanding system [39] . 
approximates the infrared regime by a zero mode that is 
coupled to a two-to-two scattering kernel. Starting from 
large initial occupancies, this description gives the scal¬ 
ing Eq. 0 in the weak coupling limit. Approximating 
the infrared dynamics by a single zero mode in a kinetic 
description can be a crude approximation. A defect of 
such an approach, known for systems in a fixed box, is 
that condensation occurs within a finite time even for 
infinite volumes [65l |66] . Since any coherent mode can¬ 
not spread out faster than the speed of light, this re¬ 
sult is an artifact of the approximation. Fortunately, the 
proper scaling of the condensation time with the volume 
can be understood in quantum field theory based on sys¬ 
tematic large-expansions HZIISH]. At next-to-leading 
order (NLO) the IR dynamics is described by an effec¬ 
tive, vertex-resummed kinetic description. This vertex- 
resummed kinetic description accurately reproduces the 
results from classical-statistical simulations [69], as has 
been shown before for non-expanding systems. 

Content of this work. We will analyze longitudinally 
expanding Wcomponent scalar field theories both nu¬ 
merically using classical-statistical simulations and ana¬ 
lytically using large-methods. We investigate the dis¬ 
tribution function /(p), whose fixed point scaling form 
is shown schematically in Fig. with different charac¬ 
teristic momentum ranges. There is a significant range 
of momenta above a characteristic screening scale, with 
large but still perturbative occupancies, where 1/A ^ 
f{p) 1. This is the regime where comparisons to stan¬ 
dard kinetic theory can apply. We note however that 
such comparisons are predicated on the interactions be¬ 
ing sufficiently local in momenta. That this is the case is 
not assured and understanding whether our results can 
be explained using kinetic theory is an important objec¬ 
tive of our work. 

In order to arrive at a deeper understanding of the dy¬ 
namics sketched in Fig. we will expand significantly 
on our previous work (Tj. We will study in detail a 
further scaling regime at hard momenta that was only 
discussed cursorily in our previous work. The scaling 
form of the distribution function exhibits an additional 
infrared regime with nonperturbative occupancies ^ 1/A. 
As noted, this regime is not well described by standard 
kinetic theory nor in prescriptions assuming an additional 



FIG. 2. Schematic scaling form of the scalar occupation num¬ 
ber distribution / as a function of (transverse) momentum. 
Classical-statistical field theory is valid for / ^ 1. Kinetic 
theory is restricted to / ^ 1/A if interactions are suffi¬ 
ciently localized in momenta. In scalar field theory, a ‘vertex- 
resummed’ kinetic description from a 1/A-expansion at NLO 
can also describe the IR momentum range. 


single zero mode. Most importantly, it turns out to have 
significant support in integrals for bulk quantities such 
as the longitudinal pressure Pl- Remarkably, for the 
scalar field theory, we will show that the effective vertex- 
resummed kinetic description from a l/A/'-expansion at 
NLO describes the physics over the entire IR momentum 
range. 


After specifying the model and parameters in Sec. |n| 
we will determine in Sec. Ill the universal exponents and 
scaling functions of the nonthermal fixed point of the 
longitudinally expanding scalar theory. We show that 
the system at low momenta behaves as a nonrelativistic 
superfluid Bose gas with an isotropic distribution func¬ 
tion but distinct scaling exponents due to expansion. In 
particular, we analyze how condensation occurs and find 
that the condensation time increases with volume. We 
further compute the properties of the scaling regimes at 
higher momenta and verify their insensitivity to varia¬ 
tions of the initial conditions and microscopic parame¬ 
ters. We discuss the impact of the different momentum 
ranges on the energy and pressure evolution of the ex¬ 
panding fluid in Sec. HYi We demonstrate that the de¬ 
viations from standard kinetic descriptions indeed arise 
from the infrared regime contributing significantly to the 
longitudinal pressure evolution for large initial occupan¬ 
cies. In Sec. |V| we estimate for the longitudinal pressure 
the nature of the quantum corrections and the time at 
which these become relevant for the evolution. A sum¬ 
mary is provided in Sec. |Vl| Some technical aspects of 
the computations are outlined in the Appendices 011 
and [cl 
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II. SCALAR FIELD THEORY WITH 
LONGITUDINAL EXPANSION 

We will consider here the nonequilibrium dynamics of 
scalar field theories with longitudinal expansion. This is 
most conveniently done in terms of Bjorken spacetime co¬ 
ordinates using the (longitudinal) proper time r, longitu¬ 
dinal rapidity r] and transverse coordinates 
The relation to Minkowski coordinates and is given 
by 

r = y/ {x^y — {x‘^Y 5 V = atanh • ( 2 ) 

The spacetime metric in Bjorken coordinates takes the 
form 

= diag(l,-l,-l,-T2), (3) 

which characterizes one dimensional expansion in the lon¬ 
gitudinal direction with the expansion rate 1 /r. 

Within this setup, we will consider a massless N- 
component scalar field theory with 0{N) symmetry and 
quartic interaction defined by the classical action, 

5'b] = J d'^X • 

(4) 

Summation over the a = 1 , 2 , ... , A scalar field compo¬ 
nents is implied here and the integration proceeds over 
f d^x = f dr d^XT dr] \/—g{x) with \/—g{x) = r result¬ 
ing from the metric determinant g{x) = det{g^y{x)). 

In the corresponding quantum theory, the fields are 
replaced by scalar Heisenberg field operators 4>a(xT, 77 , r) 
and we denote their expectation value by 

(pair) = {^a{'X.T,V,x)) ( 5 ) 

for spatially homogeneous systems. To make contact 
with kinetic descriptions, we shall consider the single par¬ 
ticle distribution in spatial Fourier space, /{pt^Pz^^)- 
This quantity may be extracted from the connected part 
of the anticommutator expectation value [42] 

Fab{xT-^T,V-V',X, t') 

= ^({^a(xT,?7,T),$6(Xr,7?',T')}) - (paMcpbix') . (6) 

The distribution function can be defined us¬ 
ing the spatial Fourier transfornj^ of the equal 


^ The Fourier transform is given by 

FabiPT, I', T, t') = j (fxTdri Fabi^T, V, T, t') + ^ 

where the rapidity wave number u is the conjugate momentum 
variable to 77 . 


time statistical correlation function F{pT,i^,r) = 
Faa{pT^ ^5 ^5 t) /a, its temporal derivatives 

F{pT,y,T) = tt' drd'^F{pT,y,T,T')\r=T' and 
F{pT,iy,r) = {rdr F r'd'^) F{pT,iy,r,r')\r=r>/2 ac¬ 
cording to [42] 

f{PT,Pz^r) + i 

= ^JF{pT,v,T)F{pT,iy,T) -F‘^{pt,v,t) . (7) 

The longitudinal momentum variable is identifiecQ as 
Pz = T^/r. 

We employ two different sets of Gaussian initial con¬ 
ditions formulated in terms of the coherent field expec¬ 
tation value 0 a (^ 0 ) and the single particle distribution 
f{PT,Pz^^o) at the starting time tq. The first set is cho¬ 
sen to exhibit a large characteristic occupancy with 

fiPT,Pz,To) = y + , (8) 

and vanishing coherent field 0a('^o) = 9r0a('ro) = 0. Here 
the parameters describe the initial occupancy no, initial 
anisotropy and the characteristic momentum scale Q 
of the distribution at initial proper time. We will refer 
to this as ‘overoccupation’ initial conditions. 

Our second set of initial conditions, are coherent field 
initial conditions, characterized by a large coherent field, 

(pai-ro) = cro ]l ^ Sal , dr(pa(jo) = 0 , (9) 

and vanishing occupancy /{pt^Pz^^o) = 0- Details of the 
numerical implementation of initial conditions are well 
documented in the literature [laiToiiii] and are summa¬ 
rized in App. 

Both sets of initial conditions admit a rigorous descrip¬ 
tion of the dynamics in terms of classical-statistical field 
theory to leading order in the coupling A dziiiniiiiiizii- 
[75]. Since we are interested in the weak coupling limit, 
with A ^ 0+ but A/ finite for typical momenta, the map¬ 
ping is valid at all simulation times. 

For the overoccupation initial conditions (|^, the 
Monte Carlo sampling and solution of the classical field 
equation leads to a well defined continuum limit of the 
corresponding quantum theory m- A rigorous mapping 
of the quantum dynamics with coherent field initial con¬ 
ditions actually involves two steps, which is some¬ 
times overlooked in recent literature. One first separates 


This choice of variable is motivated by the form of t he longitu- 
dinal kinetic term {d^(pa)l'F^ in the field equation ( |l 0 | > below. 
In App. ^ we will explain in greater detail how to compute the 
correlation function numerically and shall discuss an alternative 
definition of the single particle distribution that is also employed 
in the literature |20lI70| . 
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Transverse momentum: pj / Q 

FIG. 3. Snapshots of the occupation number distribution at 
vanishing longitudinal momentum at different times for the 
coherent field initial conditions in Eq. with coupling con¬ 
stant chosen to be A = 10“^. Starting with the pure vacuum 
(gray) distribution, fluctuations are built up via the para¬ 
metric resonance mechanism (red circles) and lead to an ove¬ 
roccupied system (blue triangles). The latter approaches a 
nonthermal fixed point (green rhombi and black squares). 

the field into a large coherent part and a small fluctu¬ 
ation and linearizes with respect to the latter. Though 
small initially, the fluctuations grow exponentially be¬ 
cause of an instability muzz]. Once they become larger 
than unity, but are still much smaller than the coherent 
field part, one uses the result as the input for a classical- 
statistical description which is fully nonlinear. The whole 
procedure has a well defined continuum limit such that 
one gets the physical quantum result im [III ITT]. In 
the following, we will employ the coherent field initial 
conditions of Eq. § for finite values A ^ 1 and mo¬ 
mentum cutoffs only to show that after the parametric 
time ^ log^/^(A“^) the initial field decays [TT] and one 
recovers the fluctuation initial conditions of Eq. ^ . 

We solve numerically the nonlinear classical field equa¬ 
tions of motiorj^ 

(^r + 7 - 4 - ‘Pa=0, (W) 

using leapfrog integration. We discretize the equations 
on anisotropic lattices of size x with lattice spac- 
ings ut, CLf]^ and use a time step a^- which is small 
compared to all other dimensional parameters, namely, 
ar ^ to guarantee the stability of the algo¬ 

rithm. We have explicitly checked that our results are 
insensitive to changes of the discretization parameters. 


^ The coupling dependence can be conveniently scaled out of the 
equations with cp —)► (pjy/X, thus entering only the initial condi¬ 
tions. 


In Fig.[^ we follow the time evolution for coherent-field 
initial conditions in Eq. (§. These initial conditions lead 
to a parametric resonance instability at early times in 
the static case [78H8Q] as well as the expanding scalar 
field theory HD [771 despite the longitudinal expansion. 
As shown in Fig. this instability results in rapid parti¬ 
cle production (red circles) and an overoccupation of low 
momentum modes built up on a short time scale (blue 
triangles) while the coherent field decays. In this case, 
we define the characteristic momentum scale as Q = do 
of Eq. and identify an overoccupied system with a 
similar form as given by Eq. (§. Hence starting with 
an overoccupied single particle distribution may be re¬ 
garded as an already time evolved system with coherent 
field initial conditions. At later times (green rhombi and 
black squares), the system approaches a nonthermal fixed 
point, whose universal properties we will analyze in the 
following. 

We note that a similar argument can be made in 
the case of the longitudinally expanding non-Abelian 
plasma. Starting from the Glasma, which consists of 
strong boost invariant color fields Hi], particles are 
rapidly produced via plasma instabilities jSTHSb] and lead 
to overoccupation [ 86 ]. An overoccupied non-Abelian 
plasma approaches the nonthermal fixed point found in 
Refs, [m in the weak coupling limit. 

III. NONTHERMAL FIXED POINTS OF 

LONGITUDINALLY EXPANDING SCALARS 

Classical-statistical field theory efficiently describes 
the far from equilibrium dynamics of highly occupied 
many-body systems. Numerical simulations in this con¬ 
text reveal that highly occupied systems across a wide 
range of energy scales spanning inflationary cosmol¬ 
ogy dll dS ED], heavy ion collisions [2TH26] . ultracold 
atoms |2D[291|3D and, as has been recently suggested, 
holographic superfluids [87] display universal behavior 
characterized by nonthermal fixed points. 

In the vicinity of a nonthermal fixed point, strongly 
correlated quantum systems lose their sensitivity to the 
initial conditions and other microscopic parameters of 
the theory EDIllESI- The single particle distribution 
function can then be divided into separate scaling regions 
in momentum space, each of which follows a self-similar 
evolution 

f{PT,Pz,T) = T°‘fs{T‘^PT,T'^Pz), (H) 

where r is the proper time. Different aspects of the dy¬ 
namics can be described in terms of the universal scaling 
exponents a, P and 7 and (proper) time independent 

scaling function pl lMlISS] . 

In general a nonthermal fixed point can exhibit multi¬ 
ple scaling regions simultaneously, with a different scaling 
behavior realized in each of these separate momentum 
regions. Such multiple scaling behavior can emerge as 
a consequence of multiple conservation laws such as the 








6 


self-similar transport of a different conserved quantity 
such as energy or particle number in one regime com¬ 
pared to another. Another possible reason includes a 
change of the underlying dynamics, for example, from 
a perturbative to a nonperturbative mechanism below a 
characteristic momentum scale. In the context of wave 
turbulence, each momentum region is usually referred to 
as an inertial range of momenta, with various examples 
provided in the literature [iniETlIMlIMlIMlESlS]. 


Within the classical-statistical simulations of the ex¬ 
panding scalar theory in [T], we identified three distinct 
inertial scaling regimes: i) An inverse particle cascade 
at very soft momenta with large occupancies / > 1 /A, 
ii) an intermediate range of momenta extending dynami¬ 
cally up to momenta close to the hard momentum scale Q 
of the system, and hi) a hard momentum regime for mo¬ 
menta Pt ^ Q- In the following, we will expand on the 


discussion of the longitudinally expanding scalar theory 
in [T] and present detailed results uncovering the proper¬ 
ties of the three inertial regimes. 

Since the nonthermal fixed point structure is insensi¬ 
tive to the details of the initial conditions, we will first 
focus on the time evolution for a single initial condition 
and then verify that the results are insensitive to the ini¬ 
tial conditions (summarized in Table further below for 
parts of the attractor. Unless stated otherwise, we con¬ 
sider a four component (0(4)) scalar theory with overoc¬ 
cupation initial conditions, as in Eq. (|^, with an initially 
isotropic system (^o = 1 ) and with the occupancy param¬ 
eter no = 35 at the initial time Qtq = 10 ^. All quantities 
in this paper are normalized by appropriate dimensional 
scales. If not written explicitly, we imply r i-^ t/tq and 
Pi ^ Pit which also holds for the scaling behavior in¬ 
troduced in Eq. (11). 

Our results for the transverse momentum dependence 
of the spectrum at different proper times in the evolution 
are compactly summarized in Eig. which shows snap¬ 
shots of the single particle distribution /{pt^Pz = 0 ^^) 
for vanishing longitudinal momenta. One clearly observes 
the emergence of three distinct scaling regimes, at small, 
intermediate and high transverse momenta. We reported 
in [ 1 ] that the low momentum behavior - indicated as 
i) - is characterized by a strong infrared enhancement 
with an approximate power law dependence while 
the intermediate momentum regime - indicated as ii) - 
features an approximate power law shown by the 
gray and black dashed lines. 

The inertial ranges of momenta for the scaling regimes 
i) and ii) are seen to shift dynamically towards lower 
momenta as a function of time as indicated by the pur¬ 
ple dashed lines in Eig. At late times, an additional 
scaling window emerges at momenta Pt ^ Q- This novel 


regime was unanticipated and we shall discuss it at length 
shortly. In this regime - indicated as iii) - the single par¬ 
ticle distribution becomes approximately independent of 
the transverse momentum up to very high transverse mo¬ 
menta where the spectrum shows a rapid fall off. 

Each of the three momentum regions corresponds to a 



Transverse momentum: pj / Q 

FIG. 4. Dependence of the single particle distribution on the 
transverse momentum (for vanishing longitudinal momentum 
Pz = 0) at various times. One observes the emergence of 
three distinct scaling regimes labeled as i), ii) and iii) and 
separated by purple dashed lines. Each regime corresponds 
to the inertial range of momenta for a nonthermal fixed point 
as discussed in detail below. The black and gray dashed lines 
illustrate the power law dependence in regimes i) and ii). 


different self-similar evolution of the single particle dis¬ 
tribution with the scaling behavior in Eq. (11). The ex¬ 
ponents determine the evolution of the typical amplitude 
/typ ^ and of typical momenta PT,typ ^ and 
Pz,typ ^ Moreover the signs of ^ and 7 provide 

the effective direction of transport of particles or energy 


Cc 

^0 

)nfig 

no 

uratic 

Qto 

>n 

N 

Nt 

Lattice pai 
Nr, 

:ameters 

Qar 


0.5 

5 

1000 

4 

128 

768 

1 

5 

1 

4.5 

1000 

4 

128 

768 

1 

5 

1 

5 

1000 

4 

128 - 256 

512 - 4096 

0.5-1 

2.5-5 

1 

6 

1000 

4 

128 

768 

1 

5 

1 

7.5 

1000 

4 

128 

768 

1 

5 

1 

15 

1000 

4 

192 

768 

0.67 

5 

1 

35 

1000 

4 

96 - 384 

384 - 768 

0.25-2.5 

4-10 

1 

80 

1000 

4 

128 

768 

0.4 

5 

2 

15 

1000 

4 

128 

768 

0.85 

5 

4 

15 

1000 

4 

128 

1024 

1 

5 

1 

35 

100 

4 

192 

1024 

0.75 

25 

1 

5 

1000 

2 

192 

768 

0.5 

5 

Eq. (91 

1000 

4 

192; 256 

256 - 768 

0.2-0.9 

6-20 


TABLE I. Discretization parameters for different initial condi¬ 
tions. Unless stated otherwise, these parameters are employed 
for all classical-statistical lattice simulations. In the follow¬ 
ing, the configurations of the last three lines are referred to 
as “Qto = 100”, “(!^(2) scalars” and “Coherent field IC” (c.f. 
Eq. §), respectively. 
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FIG. 5. Time evolution of the condensate observable with 
overoccupation initial conditions, plotted for different vol¬ 
umes V = VrLr). The lattice spacings for the smallest volume 
are Qar = 0.75 and = 5 • 10“^ on a 64^ x 384 grid. At 
sufficiently late times, the results become independent of the 
volume, signaling the formation of a condensate. The con¬ 
densation time increases with volume. For large volumes, the 
condensate growth follows power laws that are marked by 
a dashed black line. A similar late-time behavior is found 
for other initial parameters, as for an earlier initial time 
Qto = 100 . 


within the inertial range of momenta. We will extract 
the universal exponents and the scaling function of each 
scaling region in subsequent sections. 


A. Infrared cascade and Bose condensation 

We will first present a detailed analysis of the iner¬ 
tial range at low momenta, marked by i) in Fig. [^with 
an approximate power law in the transverse momen¬ 
tum spectrum. In the study of scalar field theories with 
static geometry, such a low momentum region is associ¬ 
ated with an inverse particle number cascade of an ef¬ 
fectively massive theory and leads to the formation of a 
Bose condensate [naini]. 

We can similarly investigate the presence of a dynam¬ 
ically generated Bose condensate in the longitudinally 
expanding system by studying the finite volume scaling 
of correlation functions of the scalar field. In an ap¬ 
proach similar to that considered previously for the non¬ 
expanding theory, we examine the properties of the sta¬ 
tistical correlation function Fo(t) = F{pt = 0,z^ = 0,r) 
of Eq. While the fluctuations at finite momenta are 
essential for the formation of the condensate via an in¬ 
verse particle cascade, a volume independent correlation 
function at zero momentum Fo{t)/(VtL^) = 0 o(^) sig¬ 
nals the existence of a Bose condensate. Our results for 
this condensate observable are shown in Fig. where we 
compare the time evolution for different volumes. Start¬ 
ing from overoccupation initial conditions, one observes 


a rapid rise of the zero momentum mode at early times. 
Eventually the condensate observable becomes volume 
independent thereby signaling the formation of a con¬ 
densate. As we noted in our previous letter [T], the con¬ 
densate then decreases as Fq{t) ^ 

The condensation time Tc can be viewed as the on¬ 
set time of volume independent evolution. In our case, 
the rescaled zero mode {r — tqY^^Fq becomes constant in 
time (as marked by the constant gray dashed line). While 
for smaller volumes the condensate is created during the 
simulation time, for larger volumes its formation is not 
completed. Nevertheless from the ordering of the curves 
one observes that the condensation time increases with 
volume. While conventional kinetic descriptions have a 
finite condensation time even for infinite volumes |65l | 66 ] , 
the fact that the condensation time diverges with increas¬ 
ing volume is expected from causality and is an important 
feature of our nonperturbative numerical approach. 

We note that for the longitudinally expanding system 
causality also forbids the formation of a Bose-Einstein 
condensate if the rapidity range is too large. This is be¬ 
cause at any given time r, two points which are separated 
by a larger rapidity range than Ar]c = 21 n(rc/ro) in the 
longitudinal direction are not causally connected. Never¬ 
theless one may still observe the emergence of long range 
correlations for a system with sufficiently small rapidity 
extent A 77 < Arjc and all our considerations are only valid 
in the rapidity range where we do observe condensation. 

For nonexpanding systems, both relativistic and non- 
relativistic lattice simulations also indicate an increas¬ 
ing condensation time with volume, as tc r\j m. 

This results from the condensation behavior Fq ^ a 
consequence of the inverse particle cascade from the low 
momentum region that populates the zero mode [OTfi95] . 

As can be seen in Fig.j^for larger volumes, the conden¬ 
sate observable in the expanding theory r^/^Fo(r) /V also 
follows an approximate power law behavior in time. We 
expect from analytic considerations detailed in Sec. |mB 


that the power law exponent approaches unity at late 
times. This leads to an increasing condensation time with 
volume Tc ^ E, as we will discuss in detail in Sec. IIIIB 


Indeed the results of Fig. j^seem to support this behav¬ 
ior although with sizable discrepancies at earlier times 
and for smaller volumesA similar power law evolu¬ 
tion can also be observed in simulations with different 
initial parameters and sufficiently large volume, as seen 
for the system with the earlier initial time Qtq = 100 for 
comparison. 


® Discrepancies from Ar^/^Fo(T)/F Ar at earlier times with 
Ar = r — To are expected due to finite time effects when the sys¬ 
tem is not yet close enough to the nonthermal fixed point. Such 
early time effects can also be observed in the scaling behavior of 
other quantities as shown explicitly in Fig. Similarly, simula¬ 
tions with small volumes do not show the power law behavior at 
all because the low momentum region is not completely formed 
and the system is thus too far from the attractor on the time 
scales of condensation. 









Momentum: |p|/Q 



FIG. 6. Effective dispersion relation at different times. Devi¬ 
ations from the ct; = |p| behavior show that a mass gap is gen¬ 
erated for low momentum excitations. The time dependence 
of the effective mass is shown in the inset, extracted from 
the zero momentum frequency via |p|2 fits to the dis¬ 

persion relation and by measuring the oscillation frequency 
of the zero momentum mode. The effective mass decreases 
according to a power law. 


In the scalar field theory for a static box, both the 
condensate and a dynamically generated mass were ob¬ 
served in lattice simulations [69]. We have also checked 
the existence of a dynamically generated mass for the 
longitudinally expanding case. In analogy to the static 
box case |42j |69| , we define 


^{Pt,Pz,t) 




( 12 ) 


to estimate the effective dispersion relation in the 
medium. We find that up to statistical fluctuations, the 
dispersion relation is independent of the polar angle 0 = 
/Pz) ill the entire momentum space and only 
depends on the absolute momentum |p| = \/p% 

Our results for ci;(|p|, r), averaged over the angle are 
shown in Fig.j^at different times. A comparison with the 
(dashed) lines shows that the effective dispersion relation 
is well described by a relativistic quasi-particle dispersion 


‘^iPT,Pz,T) = \/m?(T) + |p|2 , (13) 


with a time dependent effective mass m{r). 

The time dependence of m{r) is shown in the inset of 
Fig. 1^ where we compare m{r) extracted from the dis¬ 
persion relation to the one extracted from the oscillation 
frequency Co’c(t) of the unequal-time correlator at zero 
momentum F{pt = 0, z/ = 0,r, tq). Both measurements 
agree and show that the dynamically generated mass de¬ 
creases as a function of time according to 


m(r) T . 


(14) 


FIG. 7. Distribution function for different polar angles 
0 = aicta,n(pT/pz) as a function of the momentum |p| = 
\/Pt + pi at different times t/tq = 5, 10, 20. While the dis¬ 
tribution is highly anisotropic at high momenta, the infrared 
sector remains approximately isotropic. 


Because of the large occupancy at low momenta (re¬ 
gion i) in Fig. the mass is expected to be domi¬ 
nated by modes in the infrared scaling region. Using 
^ f d^pf/uj to estimate the evolution of the mass, 
one has for infrared modes 


m 


n 

m ‘ 


(15) 


Particle number conservation n ^ (which will be 
demonstrated later in this section) leads to m ^ as 

observed in (14). To arrive at Eq. ([T^ we used uj{p) 


m which is true for the infrared modes with |p| < m. 
The contribution of the condensate to the effective mass 
scales in the same way. In contrast, one can estimate 
the contribution from the harder momentum regions ii) 
and iii) in Fig. as ^ f d^pXf/uj ^ which is 

subleading because of its faster decay in time. 

We note that this is quite different from thermal equi¬ 
librium where the effective mass is dominated by hard 
momenta at the temperature scale. The difference can 
be understood as follows. Close to the nonthermal fixed 
point, the low momentum region is enhanced and its con¬ 
tribution to the mass decreases slower than at hard mo¬ 
menta, dominating at late times. In contrast, no such 
enhancement at low momenta exists in thermal equilib¬ 
rium and all contributions to the mass are stationary. 

To better understand the nature of the inertial range 
at low momenta, we analyze the dependence of the sin¬ 
gle particle distribution /{pr^Pz^^) on the polar angle 
0 = aict8iii{pT/Pz)‘ In Fig. we show the phase-space 
distribution for different valu es of 0 as a function of the 
absolute momentum |p| = While the spec¬ 

trum is highly anisotropic at high momenta, one observes 
that the strong enhancement of the infrared region is an¬ 
gle independent to within the accuracy of the numerical 
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Rescaled momentum: (t/tq)^ Pt^Q 



FIG. 8. Rescaled distribution function of the expanding scalar 
system at low momenta as a function of the rescaled trans¬ 
verse momentum for different times. Spectra at different 
times collapse onto a single curve, demonstrating the self¬ 
similarity of the infrared cascade. For comparison, the origi¬ 
nal distributions without rescaling are shown in the inset. 


simulations. Therefore typical transverse and longitudi¬ 
nal momenta evolve isotropically, which gives 

P = 1 ( 16 ) 

for the scaling exponents defined in Eq. ( pTj ), in the IR 
inertial momentum region. 

When comparing the values of the mass m(r) in Fig.|^ 
with the typical momenta for the low momentum fixed 
point in Figs. [^and[^ we find that the low momentum 
inertial range is always bounded by the effective mass 
|p| ^ m{r). Hence the dynamics of this region is gov¬ 
erned by nonrelativistic physics with the time dependent 
dynamically generated mass m(r). Moreover from Figs.|^ 
and one observes that the low momentum region en¬ 
tirely involves nonperturbatively large occupation num¬ 
bers with at least A/ > 5. 

We now turn to the universal quantities a, /3, 7 and the 
scaling function fs{PT^Pz) of the low momentum scaling 
region. To access the scaling regime earlier, we increase 
the initial amplitude to no = 125. Since we need to be 
sensitive to the low momentum region, we use a large 
lattice with 768^ x 512 lattice points and lattice spacings 
QaT = 2 and = 1.5 x 10“^. 

In Fig. we show the rescaled transverse momentum 
distribution f{pT^Pz = 0 ?^) a function of the 
rescaled transverse momentum r^pr- The inset gives 
the curves at different times without rescaling for com¬ 
parison. With a = 1 and /3 = 2/3 the rescaled curves at 
different times lie remarkably well on top of each other. 
According to Eq. ( pT| , this reveals a self-similar evolu¬ 
tion of the distribution function fs. If particle number 
is conserved, one has rn = tN f d^p/{27r)^f = const 
which requires a — 3/3 = —1 consistent with the values 
we used. The movement of the region to lower momenta, 


FIG. 9. The exponents /3 and a — 3/3 extracted at different 
reference times Tref. The gray dashed lines correspond to the 
values 2/3 and —1 extracted from a scaling analysis within a 
vertex-resummed kinetic description. 



Normalized scaling variable: pj / b 


FIG. 10. The normalized fixed point distribution fipr^Pz = 
0) at time 2tq. Also shown is the distribution function for a 
two component scalar field theory in a static box, previously 
considered in Ref. [69]. Both spectra can be well described 
by the functional form given by Eq. (18) as indicated by the 
black dashed line. 


as can be seen in the inset of Eig.|^ is consistent with the 
positive value of the exponent /3 > 0 observed and can 
be interpreted as being due to an inverse particle cascade 
feeding the Bose condensate at zero momentum. 

To gain insight into the statistical and systematic er¬ 
rors for the exponents, we follow Ref. [69] and investigate 
how the values of the exponents a and P depend on the 
reference time r^ef at which we start our self-similarity 
analysis. To this end, we compare the rescaled distribu¬ 
tion function at r^ef with the rescaled distribution func¬ 
tion at different times up to r/rref < 5, and extract the 
best values for the exponents together with an estimate 
of the error (see Refs. [22] [69] for more details on the 
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method). 

In Fig. we show the extracted values for p and the 
combination a — SP as functions of the reference time. 
Also shown are the values —1 and 2/3 as gray dashed 
lines; the former corresponds to particle number conser¬ 
vation and the latter to an elastic scattering process with 
a resummed matrix element. (This effective elastic scat¬ 
tering interpretation of the numerical data will be dis¬ 
cussed shortly in Sec. [Ill B] ) One finds that both of these 
quantities approach the gray lines. We thus have 

a = l, /3 = 7 =? (17) 

in the vicinity of the nonthermal fixed point with an error 
estimate of about 10 % for /3 and typically 10 — 15% for 
a —3/3^ as can be determined from the error bars in Fig.|^ 
We can also specify the functional form in the scaling 
region as 

XfsiPT,P.) = (|p|/ 5 )«<+(|p|/ 5 ).> > (18) 


where the non-universal constants a and b can be fixed 
by the constraints A/^dpj = b) = a/2 and A/^(|p| = 
b) = —a(/T:< + /t:>)/( 46). In Fig. 10 we show the normal¬ 
ized distribution function f{pT^Pz = 0) at r = 2ro as 
well as the corresponding distribution function for a two 
component scalar field theory in a static box, discussed 
at length in Ref. [69] . The scaling functions of both the 
static and the longitudinally expanding system are given 
by Eq. ( jlsj ). We used the values = 0.5 and =4.5 
for the scaling function in the figure but we observe that 
/^> slowly grows with time approaching the value ^ 5 
expected also for nonrelativistic systems [271 Eoj. 

Indeed Ref. [69] showed that the IR dynamics of a 
scalar field theory in a static geometry is governed by 
effectively nonrelativistic dynamics with a nonperturba- 
tive mechanism that describes an inverse particle cas¬ 
cade. While the IR scaling function in Eq. (18) can be 
observed in both expanding and nonexpanding scalar sys¬ 
tems, the temporal exponents a = 3/2 and ^ = 1/2 
for the static box in three spatial dimensions differ from 
the corresponding values oi a = 1 and (3 = 2/3 for the 
expanding theory. This difference stems from the dif¬ 
ferent geometry; in the expanding case, one has a time 
dependent effective mass and particle number conserva¬ 
tion gives rn = const instead of n = const for the fixed 
box. However these differences do not alter the underly¬ 
ing dynamics in a fundamental way. We will show that 
the IR dynamics of the expanding system, like that of 
the static case, is consistent with a description in terms 
of a vertex-resummed kinetic theory. 


B. Low momentum region from vertex-resumed 
kinetic theory 

We will now extract the exponents a, p and 7 for the 
inertial range at low momenta from a vertex-resumed ki¬ 
netic theory and show that their values are consistent 


with our numerical results. We will subsequently relate 
these universal exponents to the dynamics that leads to 
the formation of the Bose-Einstein condensate. 


1. Universal exponents 


The kinetic equation for the distribution function / in 
the longitudinally expanding background reads as 

(dr - f{pT,Pz,T) = C[f]{pT,Pz,T) , (19) 

with the collision integral C[f]{pT,Pz^^) representing the 
scattering process and with the redshift term —p^ /t 
that follows from the boost invariance of the system in 
the longitudinal direction [2ll|47l|96]. 

Since the occupation numbers are nonperturbatively 
large in the low momentum region / > 1 /A, a perturba¬ 
tive approach is problematic. However for an N compo¬ 
nent scalar theory, one can employ a l/A/" expansion of 
the two particle irreducible (2PI) effective action [42ll68) . 
To NLO in this l/N expansion, the temporal dynamics 
of the distribution function can be expressed as a kinetic 
equation with the vertex-resummed collision kernel intro¬ 
duced in Refs. [69l[97]. Specifically, in this approach, all 
corrections to the coupling constant of the Wcomponent 
scalar field theory are included in an effective coupling 
Aeff. The collision integral for this effective kinetic de¬ 
scription in the infrared describes an elastic 2 2 scat¬ 

tering process with Aeff[/], 

C[f]{pi) = I ^ [(/i + f2)hh - /l/2(/3 + /4)] . 

( 20 ) 


Here we used fi = f{pT,iiPz,ii^) ^ 1 with pi = p and 
suppressed the arguments of the functionals J 
and Ag^[/]. Their precise forms are written in Ap¬ 
pendix [^ 

To extract the universal scaling exponents of the low 
momentum region, we follow Refs. pniisn] and study the 
temporal scaling behavior of each side of Eq. ( [l^ . Using 
the scaling ansatz in Eq. ( [TT| , the left hand side becomes 

(a + /3 pt dp-^ + (7 - 1) p^ dp -^) fs {pt,Pz) , (21) 

where we have abbreviated pr = t^pt and Pz = r^Pz- 
Similarly one can determine the scaling property of the 
collision integral 


C[f]{pT,Pz,r) =T^ C[fs]{pT,Pz) 


( 22 ) 


where C[fs\{pTiPz) is time independent. The exponent 
p = /i((a, /3, 7 , a) characterizes the scaling behavior of the 
collision integral, where we also included the scaling ex¬ 
ponent a of the effective mass m{T) ^ . Substituting 

Eqs. ( 21 ) and ( 22 ) in Eq. ([T^ leads to the time indepen¬ 


dent fixed point equation for fs{PT^Pz) 


oifs{PT,Pz) Pprdp^ fs{PT,Pz) 

+ (7 - dp. fs{PT,Pz) = C[fs]{pT,Pz) , (23) 
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where we have simply relabeled pt Pt and Pz ^ Pz^ 
and to the scaling relation for the temporal exponents 


and therefore integrating the kinetic equation (19) leads 
to 


a — 1 = p. 


(24) 


A scaling solution of the kinetic equation has to satisfy 
both equations simultaneously. A detailed discussion of 
the fixed point equation (23) and its solution is beyond 


the scope of this paper. We will focus on the scaling 
relation Eq. (24) since we are primarily interested in the 
scaling exponents. 

To determine the scaling exponent of the collision in¬ 
tegral /i, we first observe that since the inertial range at 
low momenta is bounded by the effective mass m(r), the 
dispersion relation becomes 


^{Pt,Pz,t) ~ m{T) + 


2m{r) 


(25) 


Because the leading mass term cancels in intermediate 
steps, one has to determine the scaling property of the 
momentum |pp. As our numerical results in Fig. ([^ 
show, the momentum distribution in the IR is isotropic. 
Hence we can assume that longitudinal and transverse 
momenta as well as the modulus all scale equally with 
^ = 7- 

The precise form of the exponent p is derived in Ap¬ 
pendix and reads as 

p = a — 2l3 -\- a ^ (26) 


where the reader will recall that the scaling exponent 
a controls the temporal evolution of the effective mass. 
From the scaling relation in Eq. (24), and using the 
isotropy condition, one obtains 


/? = 7 = 


1 + cr 


(27) 


Since our numerical results show that the mass scales 
with a = 1/3, the scaling relation deduced from the 
effective kinetic theory results in momentum exponents 
that decrease with /3 = 7 = 2/3, which is consistent with 
the numerical result we obtained from classical-statistical 
simulations. 

To determine the scaling exponent a from the vertex- 
resummed kinetic theory, we need a second condition. 
The form of the elastic scattering kernel in Eq. (19) sug¬ 


gests that particle number (defined as the integral of / 
over phase space) is conserved. One can confirm that 
the collision integral ( 20 ) vanishes when we integrat^ 
N f (PpT diy/{27r)^C[f\ = 0. Further, the corresponding 
integral over the left hand side of Eq. (19) can be iden¬ 


tified with the total time derivative of the distribution 
function 

(dr - f{PT,Pz,T) = ^f(pT,v/T,T) , (28) 


^ Recall that u = rpz- 


0 = 


ip 


(fpT dv d 

- —f{PT,v/T,T) = — Tn. (29) 


(27r)2 


dr 


Thus rn = const in the IR inertial range; from the self¬ 
similar form of Eq. (11), one then obtains 


a = 2(3 + p — 1. 


With Eq. (27), this leads to 

(a = (l-l-3cr)/2. 


(30) 


(31) 


Again, plugging in a = 1/3 we obtain a = 1 , which 
is again consistent with the result obtained from the 
classical-statistical numerical simulations. One may 
therefore conclude that the vertex-resummed kinetic the¬ 
ory obtained from the NLO 1/N expansion in the IR 
inertial range is consistent with first principles classical- 
statistical numerical simulations of the expanding scalar 
theory. 


2. Condensate formation 


With the scaling properties of the IR inertial region, 
we are now able to further analyze the power law growth 
of the zero mode Tb(r) = F{pt = 0,p^ = 0,r) in Fig. 
for sufficiently large volumes. The correlation function F 
can be decomposed into a particle and condensate part 
according to 


F{pt,v,t) 


f{PT,Pz,T) 

T oj{pt,Pz,t) 


+ {2TrfS^^'>{pT)S{u)(j)l{T). 

(32) 


For finite momenta, the “particle” part (c.f. App. H 
of this expression provides an alternative definition of 
the distribution function. Our expression in Eq. ( [^ is 
essentially the same as the one derived previously for a 
static box [69], except for the additional factor 1/r due 
to the expanding metric. We have checked numerically 
that it converges to the previous definition in Eq. 0 with 
increasing accuracy in time. 

When considering overoccupation initial conditions ^ 
there is no condensate initially but the inverse particle 
cascade populates the zero mode according to the relation 


n /?/3 ^ , (33) 

Uj{p = 0) 


where we have used the scaling exponent of the mass 
a = 1/3. Denoting the value of the rescaled zero momen¬ 
tum correlator {r/F q{t)/V at the initial time of 
the self-similar regime as f{p = 0,ri)/{rim{Ti)V) and 
its final value at the condensation time Tc as 0 o(^c) = 
(^c/^i)^^^ 0 o(^c)) the condensation time Tc then follows 
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from the power law dependence of the zero mode in 
Eq. as 


( <^o(Tc)nm(Tj) \ ' 

y f{p = 0,ri) j 


Hence, we have Tc ^ ^ which is the same relation 

as for the non-expanding theory [69]. Most importantly 
it diverges with increasing volume as Tc ^ E for the ex¬ 
panding theory close to its attractor. 

We can also estimate the number of particles in the 
condensate Uc beyond the onset of condensation. By use 
of the decomposition in Eq. ( j^ one finds 

ric ~ uj{p = 0) (/)o ~ , (35) 

where in the last step we used the scaling properties of 
the mass and the condensate. Hence once a condensate is 
created, its particle number is approximately conserved 
and its evolution effectively decouples from the rest of the 
system. This is certainly not true for times prior to the 
formation of the condensate, where the low momentum 
region plays an essential role. 


C. Intermediate momenta 

We now turn to a more detailed analysis of the inter¬ 
mediate momentum sector characterized by the ^ IjpT 
power law in Eig. In contrast to the soft sector, the 
spectrum at intermediate momenta is anisotropic; it is 
more efficiently described by considering the dependence 
on the longitudinal and transverse momenta separately. 

Regarding the transverse momentum dependence, we 
discussed in [ 1 ] that in this case, there is effectively nei¬ 
ther a particle nor an energy flux in transverse momen¬ 
tum. Instead the energy and particle number per trans¬ 
verse momentum mode rdn/dpT and r dejdpT are time 
independent and thus conserved at each transverse mo¬ 
mentum separately. Additionally, rdn/dpT is uniformly 
distributed over transverse momenta, 

. dn 

T A -— rsj r pt 
dpT 

up to a rapid fall-off for > 2 Q. 

The local conservation of both energy and particle 
number is a general property of the anisotropic scal¬ 
ing regions of the longitudinally expanding scalar field 
theory and has also been observed in Ref. [ 21 ] for the 
non-Abelian gauge theory. In contrast, there is no sin¬ 
gle scaling solution conserving both energy and particle 
number for isotropic systems [20]. Instead only one con¬ 
servation law constrains the scaling solution. This has 
been observed not only for the scalar theory in a static 
geometry uninsiHi] but also within the low momentum 
inertial range of the expanding scalar theory where local 
particle number conservation leads to an inverse particle 
cascade. 


/ 


^ >^f{PT,Pz, t) = const , (36) 



10 100 



Time: t/tq 


FIG. 11. Time evolution of the longitudinal momentum scale 
(cr^)(T) for different initial conditions. Top figure: After a 
short transient behavior, the dynamics becomes insensitive 
to different values of no and all curves follow an approx¬ 
imate power law dependence indicated by the black 

dashed line. Comparison to the gray dashed free streaming 
curve shows the significance of momentum broaden¬ 
ing. Bottom figure: Three of the curves show the insensitivity 
of the late time behavior to variations of the initial momentum 
anisotropy ^o- The curve labeled Qtq — 100 corresponds to a 
different initialization time and demonstrates insensitivity of 
late time results to the starting time. The plot labeled 0(2) 
scalars shows the results for two component scalars. Since 
all the other curves are for iV = 4 scalars, this curve demon¬ 
strates that late time results are independent of the number 
of degrees of freedom. 


In the language of the self-similar evolution in Eq. (11), 
the local conservation of energy and particle number in 
the anisotropic scaling regime gives rise to the scaling 
relations 


(a-7-hl = 0, /3 = 0. (37) 

With these constraints taken into account, the relevant 
dynamics in this regime is that of longitudinal momen¬ 
tum broadening. 

Momentum broadening in the longitudinal direction 
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FIG. 12. Dependence of the single particle distribution 
function on the longitudinal momentum variable v — rpz at 
Pt = Q/2 for different times of the evolution. Momentum 
broadening leads to a redistribution of the spectrum from 
lower to higher u. 


can be efficiently analyzed in terms of the characteristic 
longitudinal momentum scale 


^liPT,T) 


I dpz plf{PT,Pz,T) 
J dpzf{PT,Pz,T) 


(38) 


In Fig.we show the time evolution of {(j‘ 1 ){t) (that is 
cf‘1{pt over transverse momenta within the 
intermediate inertial range ii) for various simulation pa¬ 
rameters. The top panel shows the behavior under vari¬ 
ations of the initial overoccupancy no, while the bottom 
panel corresponds to variations of the initial anisotropy 
^ 0 , the initial time Q tq as well as the number of field com¬ 
ponents. While for different initial conditions the early 
time dynamics ranges from free streaming to an approx¬ 
imately constant behavior, one clearly observes that the 
late time behavior becomes insensitive to all of these as¬ 
pects. Ultimately all curves follow a universal scaling 
behavior in time. 

According to Eq. ( pTj ) the scaling behavior of cFz is 
characterized in terms of the scaling exponent 7 , 


<^1{Pt,t)'^t . (39) 


We have checked that this scaling holds for all transverse 
momenta within the inertial range. By performing a com¬ 
bined analysis of all data sets, we obtain an estimate of 
the scaling exponent 


27 = 0.66 ± 0.05 . (40) 

A comparison with the data in Fig. El shows that the 
observed scaling behavior is very well reproduced by an 
approximate power law dependence. 

We can also analyze the dynamics of momentum 
broadening on a more microscopic level by directly con¬ 
sidering the single particle distribution at a typical trans¬ 
verse momentum pl^^' within the intermediate inertial 


FIG. 13. Normalized distribution fs{pzl(yz) for the fixed 
point at typical intermediate transverse momenta The 

rescaled spectra of Fig.[^for no = 35 and ^0 = 1 at different 
times collapse onto a single curve, proving the self-similarity 
of the evolution (right labels). As shown, the functional form 
of the fixed point distribution is independent of the initial 
conditions (left labels) and can be described by a normalized 
Gaussian distribution (black dashed curve). 


range ii). Our results for the time evolution of the longi¬ 
tudinal spectrum are presented in Fig. which shows 
snapshots of the distribution at = Q/2 for differ¬ 
ent times. The spectrum is shown as a function of the 
rapidity wave number u = rpz on the horizontal axis, 
which effectively amounts to undoing the redshift of lon¬ 
gitudinal momenta. While for a noninteracting system 
the distribution becomes time independent when plotted 
as a function of this is clearly not the case for the 
longitudinally expanding scalar theory. Instead, Fig. \T2\ 
shows a clear broadening of the longitudinal momentum 
distribution as a function of time along with a decreasing 
amplitude /o(t) = f{pT = Pt^',P z = 0, r). 

Once we take into account the decrease of the overall 
amplitude of the distribution according to 


/o(r) -- r° 


(41) 


with a —2/3 determined by the scaling relation in 
Eq. (37) and normalize the longitudinal momenta by 


the typical momentum scale Pz /the rescaled 
distribution approaches a stationary (time independent) 
fixed point fsiPz/o'z)- This is shown in Eig. 
where we present the rescaled distribution f{pT = 
Pt^'^P z^^)/fo{^) as a function of the rescaled longitu¬ 
dinal momentum Pz/^z{Pt ^'The fact that the re¬ 
sults at different times are indistuinguishable confirms 
that the dynamics of momentum broadening can be ac¬ 
curately described by a self-similar scaling behavior as in 
Eq. 

The functional form of the fixed point distribution is 
well described by a Gaussian distribution of width az^ 
which is shown as a black dashed line. One can also 
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verify explicitly from Fig. that the functional form of 
the fixed point distribution is independent of the initial 
conditions as well as the number of field components. We 
conclude that the universal properties of the intermediate 
momentum regime can be summarized in terms of the 
scaling exponents 

a - -2/3, - 0 , 7-1/3 (42) 

and the scaling function 

- ^exp , (43) 

with a nonuniversal constant A. 

Remarkably, these features of the intermediate mo¬ 
mentum regime of the expanding scalar theory are iden¬ 
tical to those that we demonstrated previously for ex¬ 
panding overoccupied non-Abelian plasmas [2TI |22] . In 
[1] , we argued that this agreement provides strong proof 
of universality, since the underlying elementary micro¬ 
scopic dynamics of the two theories are very different. 
Our results here significantly strengthen the arguments 
presented in [T]. 


D. Hard momentum region 


As noted previously, the inertial range of momenta for 
the intermediate momentum regime shifts towards lower 
momenta. At late times, an additional scaling regime 
emerges for hard transverse momentum modes, where the 
Pt spectrum in Fig. becomes approximately indepen¬ 
dent of the transverse momentum. Since the spectrum in 
this regime is highly anisotropic (with the typical longitu¬ 
dinal momenta much smaller than the typical transverse 
momenta), we will again analyze the transverse and lon¬ 
gitudinal momentum dependence separately. Just as in 
the case of the intermediate momentum range, we find 
that there is effectively neither a particle nor an energy 
flux in the transverse direction. Thus as was the case 
in our prior discussion for the intermediate momentum 
range, the relevant dynamics is that of longitudinal mo¬ 
mentum broadening and the scaling relations in Eq. (37) 
also hold for the hard momentum sector. 

The momentum broadening in this novel hard sector 
can be efficiently described in terms of the hard longitu¬ 
dinal momentum scale 


Ai(r) = 


{{dr,d^,ipa{T, Xt,v)) {dr,d^<Pa{T,XT,ri))) 
( (d^iPair, XT, ri)) XT, ri))) 


(44) 


with summation over spatial indices fi = x^,x^,ri implied 
and where ( ) includes a volume average. 

Within the quasiparticle picture (see Eq. (32)) and us¬ 
ing u){pt,Pz,t) — IpI for hard momenta, the hard scale 
Ai can be expressed as 


_ j (fpT j dp^ pI L0{pt,Pz,t) f{pT,Pz,T) 

~ J d'^PT J dpz I^{Pt,Pz,t) f{pT,Pz,T) 




FIG. 14. Upper panel: Time evolution of the longitudinal mo¬ 
mentum scale at hard transverse momenta for different initial 
conditions. At late times all curves approach a common scal¬ 
ing behavior (black line). Lower panel: The scaling exponent 
27 is extracted from the logarithmic derivative in the lower 
panel. The gray lines indicate the value 27 = 1.05 ± 0.1 ob¬ 
tained from an average over all shown simulations. 


which is closely related to a transverse momentum inte¬ 
gral of cr^(pT)- However a crucial difference is the ad¬ 
ditional weighting by the mode energy cj, which ensures 
that is dominated by high momentum modes. 

The time evolution of the longitudinal hard scale is 
shown in the upper plot of Fig. One observes that 
after an initial transient behavior the results for different 
initial conditions approach a common power law decay. 
By use of Eqns. 0 and ( [T7| ), this scaling behavior is 
described by A| ^ 7--27 g^j^d the dynamical exponent 
7 can be extracted from the hard scale by taking the 
logarithmic derivative 


2j{t) 


dlogr 


(46) 


The resulting exponent is shown in the lower plot 
of Eig.[^as a function of time for different initial condi¬ 
tions. One observes that all curves approach a constant 
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FIG. 15. Normalized fixed point distribution fsiPzjcFz) for 
the fixed point at hard transverse momenta. The rescaled 
spectra at different times collapse onto a single curve, proving 
self-similarity of the evolution (right labels). The functional 
form of the fixed point distribution is independent of the ini¬ 
tial conditions (left labels) and can be described by a normal¬ 
ized hyperbolic secant function sech(x7r/2) = l/cosh(a:7r/2). 


value, 


27 = 1.05 ± 0.1 , 


( 47 ) 


independent of the initial condition. By use of the scaling 
relations (37) we infer from this analysis the approximate 


values of the scaling exponents 

(a-- 1 / 2 , /3 = 0 , 


7 ~ 1/2 


( 48 ) 


We note that the larger value of 7 compared to the inter¬ 
mediate momentum sector characterizes the fact that the 
momentum broadening is less efficient for high momen¬ 
tum modes. Nevertheless the observed behavior is dis¬ 
tinctly different from free streaming and significant mo¬ 
mentum broadening occurs throughout the entire time. 
Moreover we can understand the motion of the bound¬ 
ary between the intermediate and hard sectors (dashed 
purple line in Fig. [3 by comparing the decrease of the 
amplitude fo{pT^r) = f{pT^Pz = O,'?") in both regions. 
The larger value of a in the hard region relative to the 
intermediate one implies a slower decrease of its am¬ 
plitude. Therefore the boundary between both sectors 
moves to lower transverse momenta and increases the in¬ 
ertial range iii) of hard transverse momenta. 

We have also extracted the scaling function fs{Pz/<^z) 
following the same procedure outlined in the previ¬ 
ous section. Our results are summarized in Fig. 
where we show the normalized distribution function 
f{PT^Pz^^)/fo{PTj as a function of the normalized lon¬ 
gitudinal momentum Pz/(Jz{pt:^)- Our results at dif¬ 
ferent times agree with each other and demonstrate the 
emergence of self-similarity with a time independent scal¬ 
ing function. We have also included the results for differ¬ 
ent hard transverse momenta Pr > which show that 


the scaling function has the same form over the entire 
inertial range of momenta. 

One observes from Fig. that the functional form of 
the scaling function is quite different from the Gaussian 
shape observed for the intermediate momentum fixed 
point. Compared to a Gaussian, the distribution in 
Fig. features a smaller width and larger tails. This 
functional form is well described by a hyperbolic secant 
sech(x ' 7 r/ 2 ), which is also shown in the figure. 


IV. IMPLICATIONS OF SELF-SIMILARITY 
FOR BULK ANISOTROPY 


Since the scalar field theory shows a rich dynamical 
structure, with different scaling behaviors observed in 
different momentum regimes, it is interesting to inves¬ 
tigate how the interplay of the different sectors affects 
the evolution of bulk observables. We will focus on the 
bulk anisotropy of the system and investigate the time 
evolution of the ratio of the longitudinal and transverse 
pressures. 

The energy momentum tensor of the scalar field theory 
is defined as 


Tf,a = {{d^^a){da^a) “ • ( 49 ) 


The transverse and longitudinal pressures are given by 


Pt 


f Til + T22 

V 2 


Pl 


1 

7^ 


/ ^ 


( 50 ) 


and the energy density is defined to be 5 = ^ T^rIV. 

Here we have used the abbreviations ^ = J (Pxt dp 
and V = Vt Lrj. 

Before we turn to a more detailed discussion of our 
results for the ratio PlIPt presented in Figs. and 
it is useful to first consider the quasiparticle expression 
for the energy momentum tensor 
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CppT dpz Pf,Pa f{pT,Pz,T) 
(27r)3 uj{pt,Pz,t) 


( 51 ) 


with ij.,a = x^,x‘^,ri and with the rapidity wave number 
= u = rpz. We note that for a proper treatment of 
the nonperturbative dynamics at low momenta, further 
contributions should be taken into account. While the 
energy density and the transverse pressure are generally 
dominated by the hard excitations of the system, the 
longitudinal pressure in particular can receive significant 
contributions from the soft sector. We have analyzed this 
behavior by performing a comparison of the integrands 

ifiPr Pt Pt f{T,PT,Pz) 

dpTdpz (27r)2 2uj{pt,Pz,t) 

d^PL Pt pI f{'r,PT,Pz) 

dpTdpz (27r)2 !x{pt,Pz,t) 

Our results are shown in Fig. [l^ where we show a con- 
tour plot of both quantities. While the integrand of the 
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FIG. 16. Ratio of the longitudinal to transverse pressure for 
initial conditions with small initial occupancies no. The ra¬ 
tio Pl/Tt exhibits an approximate scaling as expected 

from kinetic theory. In contrast, for the larger occupancies 
shown in Fig. a much slower decay of Pl/Pt is observed 
just as in the case of non-Abelian plasmas. 


transverse pressure is peaked around transverse momenta 
of the order of the hard scale, the longitudinal pressure 
receives its dominant contribution from modes with much 
softer transverse momenta. This points to the possibil¬ 
ity that the quantities Pl and Pt can be dominated by 
modes from different inertial ranges of momenta. 


The time evolution of the pressure ratio PlIPt is 
shown in Figs. (and earlier in Fig. as a function 
of time. While in all cases the bulk anisotropy of the 
system increases as a function of time, the quantitative 
behavior on the time scales observed can be different de¬ 
pending on the initial conditions. For low enough initial 
amplitudes uq < 6 shown in Fig. 16, one observes an ap¬ 
proximate PlIPt ^ scaling. In this case, a strong 

infrared enhancement of the spectrum - corresponding 
to the scaling regime i) - is not fully developed during 
our simulations. Instead the observed behavior is consis¬ 
tent with the expectation that the longitudinal pressure 
is dominated by modes in the intermediate transverse 
momentum regime ii). 


In contrast, for the larger values of no shown in 
Fig. the strong enhancement at soft momenta (given 
by Eq. @) is fully developed and has considerable con¬ 
tributions to the longitudinal pressure. In this case, the 
combined contributions of soft and intermediate modes 
lead to a slower decrease of PlIPt in time. Interestingly, 
the approximate scaling behavior observed in this case 
is quite similar to results obtained for the non-Abelian 
gauge theory also shown in Fig. This result is con¬ 
sistent with our earlier observation for the gauge the¬ 
ory case that modes with relatively small transverse mo¬ 
menta provide a sizeable contribution to the longitudinal 
pressure (c.f. Ref. [22]). 
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FIG. 17. Gontour plot of the contributions (see Eq. ( [^ ) 
to the transverse (top) and longitudinal pressure (bottom) 
for initial occupancy parameter no = 35 at time t/tq = 50. 
While the transverse pressure is dominated by hard excita¬ 
tions, one observes that the longitudinal pressure receives its 
dominant contributions from intermediate and soft momenta. 


V. RANGE OF VALIDITY OF THE 
CLASSICAL-STATISTICAL APPROXIMATION 


In this section, we will examine the range of validity of 
the classical-statistical approximation employed through¬ 
out this paper to study the far-from-equilibrium dynam¬ 
ics of a longitudinally expanding scalar theory. Gener¬ 
ally the classical-statistical approximation can be justi¬ 
fied when the typical occupancies / of the system are 
large and quantum effects play a negligible role. In par¬ 
ticular when the initial occupancies are / ^ 1/A with 
A ^ 1, as considered in this paper, quantum effects are 
sub dominant for times rguant ^ Estimating the 

exponent q requires the evaluation when quantum pro¬ 
cesses, that are initially suppressed, become important 
with the increasing dilution of the system. 


On the level of kinetic equations, the validity of the 
classical approximation can be examined by comparing 
the classical and quantum contributions to the collision 
integral. One finds that for a perturbative 2 2 scat¬ 

tering process in scalar field theories, the classical part 
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of the collision integra]^ takes the form pTl UTl H2] 
Cciassicmbl) = I 

X ((/l + /2)/3/4 - /l/2(/3 + /4)) , (53) 

with the corresponding quantum terms 

CQuant[/](Pl) = ^ j “ /1/2) • 

(54) 

While the final state Bose enhancement factors are con¬ 
sidered to be large in the classical field description, this 
is not the case when genuine quantum corrections are 
relevantly One may therefore expect that as long as the 
classicality condition 

/(Pr,Pz,T)>l (55) 


is satisfied for typical momentum modes, the quantum 
corrections in Eq. (54) are suppressed and classical scat¬ 
tering dominates. 

However as recently argued in Ref. [39] , there is poten¬ 
tially a caveat to this argument when the system is highly 
anisotropic. This is because the classical-statistical ap¬ 
proximation restricts the phase space available for each 
individual scattering to states which are already occupied 
while the quantum terms allow one to access previously 
unoccupied states by a single scattering. Assuming that 
the typical longitudinal momenta Pz ^ of occupied 
states are much smaller than the typical transverse mo¬ 
menta pt ^ At by a factor of J ^ Kz/Kt <C 1, one 
finds that a classical scattering can only lead to longi¬ 
tudinal momenta on the order of a few times ^ A;^. In 
contrast, the quantum contribution, while suppressed by 
occupancy factors, has a larger phase space and allows for 
states with longitudinal momenta up to At ^ A^ to be 
populated by a single scattering. These could in principle 
give a large contribution to the longitudinal pressure. 

We will show here that this is not the case in the para¬ 
metric domain of validity of the classical regime. We will 
begin by first estimating parametrically the importance 
of ‘quantum’ large angle 2^2 scatterings into previ¬ 
ously unoccupied phase space - in particular for modes 
with longitudinal momenta Pz ^ At ^ A^. We will 
subsequently investigate whether the dynamics of this 
process can lead to an earlier breakdown of the classical- 
statistical approximation than suggested by Eq. (55) by 
estimating the contribution of this process to the longi¬ 
tudinal pressure. 


® At leading order in 1/N expansion, the prefactor becomes 
A^/GA/”, which is analogous to the prefactor of the vertex- 
resummed kinetic equation in (pOl. 

® Similar conclusions can be reached also for higher order n ^ m 
processes and in functional nonequilibrium approaches beyond 
kinetic theory m- 


A. Parametric estimate of quantum large angle 
scattering 


We begin by estimating the number density of particles 
at large angles, which can be defined in terms of 


Aparge—angle ) 


i 


\Pz\:>^z 


d^P f{PT,Pz,T) 


(56) 


where the integral over pz is limited to modes with longi¬ 
tudinal momenta Pz parametrically on the order of Pt 
- much bigger than the typical longitudinal momenta 
Az. We can estimate the change in time of this quan¬ 
tity by considering the Boltzmann equation for 2 2 

scatterings, which, including both classical and quantum 
contributions, takes the form 


9r-/VLarge-angle('r) ~ / d^Pi / (Pp2 / d^Ps / cfipi 

^ \pi,z\^^z J J J 

-6'^‘^\pi +P2-P3 -Pi) X 


A2 


2uj\2 uj2 2c(J3 2uJ/^ 

(1 + /l)(l + /2)/3/4 — /l/2(l + /3)(1 + Ia) 


(57) 


with the abbreviations fi = f{Pi,T^Pi,z^) ^.nd similar for 
uji. Since we require pi^z to be much larger than the typi¬ 
cal longitudinal momentum of all highly occupied states, 
he. pi^z ^ At ^ Az , longitudinal momentum conserva¬ 
tion dictates that also another longitudinal momentum is 
large, At Az , meaning that both particles fall 

into the phase space with /i, /2 ^ 1. We can thus neglect 
the Bose stimulation factors in the gain term (namely, 
(l + /l)(l + /2) r\j 1) such that only quantum processes 
contribute to occupancy at large angles. Similarly, we 
can neglect the loss term altogether since / 1/2 <C 1, such 
that our estimate becomes 

9rA^Large-angle(T) ~ / dpi^^ / d'^pi^T / d^P2 X 

\Pl,z\^^z J ^ 

(58) 


We expect the dominant phase space for this scattering 
to be when the energy of all particles is on the order 
of the hard transverse scale At; therefore parametrically 
the energy denominator uji ^ At- One thus obtains the 
phase space integral for this scattering to be 

5 TiVLarge-angle('r) ^ A [ dpi^^ [ d^P 2 ,T X 

At J\pi,z\'^\z J 

J d^P2 J d^P3 J (fp4 +P2-P3 -Pi) fifi- 

(59) 

While at first sight it may seem that this phase space 
is strongly enhanced for the large angle process, a sub¬ 
tle point is that longitudinal momentum conservation re¬ 
quires pi^z + P 2 ,z = P?,,z + Pa.z to be on the order of 
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^ which is the typical longitudinal momentum of 
highly occupied particles. Implementing this momentum 
conservation constraint leads to 


large angle particles are on the order of the transverse 
hard scale At, the contribution to the pressure takes the 
form 


f 

AAarge—angle^ J d {j)\ P 2 ) 

J d^Ps J d^Pi < 5 (wi +LO2-CO3- W4) fsfi ■ ( 60 ) 


Since the typical energy of all participating scatterers is 
of the order of At and there are no further constraints 
beyond energy conservation, the phase space for the scat¬ 
tering should be of order A‘^. This leads us to the final 
estimate 




Large 


-angle (^) ^ ^^^hard(^) • 


(61) 


where A^hard('?") denotes the typical density of hard par¬ 
ticles 

-/Vhard('r) ~ j d^p fiPz,PT,T) . (62) 


Therefore the total number of particles found at large 
angles should be of the order of the scattering rate in 
Eq. mil t imes a time scale on the order of the age of the 
systenj^q 


^T™t(r) ^ A^Large-angle(^)AT ^ A^ ^ At T. 

(64) 

which should be compared to the (perturbative) classical 
contribution 


pClassic _ (g5) 

At 

Since the total number of hard particles is approximately 
conserved during the evolution, their density decreases as 
^ 1/r due to the longitudinal expansion, such that 

-/Vhard('r) ~ ) -/Vhard(To) ~ T 

and the ratio of the quantum contribution to the pressure 
over the classical contribution is given by 

pQuant A 2 

dT_A ^ (67) 

pClassic • VO/; 


Large—angleij ~) 


\2A^hardi^) A 


(63) 


B. ‘Quantum’ contributions to observables 

With our parametric estimate of the density of large 
angle particles in hand, we will now estimate the contri¬ 
bution of such quantum effects to the longitudinal pres- 
sur €3 Since the typical longitudinal momenta of the 


An alternative way to arrive at the same parametric estimate is - 
following Ref. m - to consider the rate of large angle scatterings 
that a single hard particle undergoes in the medium. Since the 
process is not Bose enhanced, this large angle scattering rate is 
simply given by the density of scattering centers N^ard times the 
cross section 


^/yColl 

Large —angle 

dr 


' ^Large —an gle-^hard • 


The cross section for a large angle scattering of hard scalar parti¬ 
cles is parametrically crLarge-angle 
of particles at large angles therefore becomes 


A^/A|.. The total number 


^^yColl .2 

-'VLarge-angle('r) ~ iVhard(T)-^ _Jv2^^j(t)t 

which agrees with the more detailed estimate in the text. 

In general, the relative importance of the leading classical and 
subleading quantum contributions (as well as other higher or¬ 
der effects) depends on the observable at hand and needs to be 
evaluated for each observable. While certain ultraviolet sensitive 
observables may receive large corrections from higher order pro¬ 
cesses, the dynamics of the bulk may be perfectly well described 
by the leading order classical process. 


hence 


At early time^^ r ^ tq, A^ ^ A|., 

pQuant^pCiassic ^ which is Very small for A ^ 1. 
With increasing time though the A^ will decrease rela¬ 
tive to as illustrated in Fig. 18 Since on large time 
scales the longitudinal pressure is dominated by modes 
in the intermediate rather than hard momentum region, 
we will assume for this estimate a la BMSS that 


Kz ^ . 


( 68 ) 


One then finds that the quantum contribution to the 
pressure becomes of the same order of magnitude as the 
classical contribution on the time scale 

AT^Quant ^ A“^/^ . (69) 

Interestingly however since /{pr^Pz ^ Az) ^ (Att)“^/^ 
in this intermediate transverse momentum regime - the 
time scale obtained from Eq. ( [6^ is parametrically the 
same as that when the occupancies in this regime be¬ 
come of order unity. We may therefore conclude that the 
‘large angle’ quantum contribution to the pressure be¬ 
comes of the order of the classical contribution only on 


Even if the initial conditions are such that the typical ~ 0 
at the instant of the collision, modes at the hard scale are very 
quickly populated by the instability resulting from the paramet¬ 
ric resonance. As noted previously, the characteristic time scale 
for this is r = To ln^/^(A“^). How rapidly this occurs is clearly 
visible from Fig. A similar argument applies for gauge theo¬ 
ries [22] [86]. 
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FIG. 18. Shown are the small angle and large-angle contri¬ 
butions to the longitudinal pressure as described in the text. 
While the large-angle contribution is suppressed initially by 
1/A, it becomes important when ~ AA|-, which occurs 
after the time . 


the time scale when the classical-statistical approxima¬ 
tion is breaking down anyway. This is indeed therefore a 
self-consistent criterion for when classical dynamics can 
be trusted. Beyond the time scale rguant, the evolution of 
the system is expected to change qualitatively. While for 
non-Abelian gauge theories there has been a lot of recent 
progress in understanding the dynamics of the quantum 
regime in more detail [38j [51] , a similar analysis for the 
scalar theory is presently not available - we will return 
to this point in the future. 


VI. SUMMARY 

In this paper, we reported on the results of classical- 
statistical simulations of the far from equilibrium dynam¬ 
ics of a highly occupied longitudinally expanding 0{N) 
scalar theory. This theory features remarkably rich dy¬ 
namics which is described simultaneously by three dis¬ 
tinct inertial scaling regimes i) at soft, ii) intermediate 
and iii) at hard transverse momenta. The soft momen¬ 
tum region is governed by an inverse particle cascade 
feeding a Bose condensate while the intermediate and 
hard momentum regions describe longitudinal energy and 
particle transport. Within each inertial range of mo¬ 
menta, the single particle distribution exhibits a self¬ 
similar behavior characterized by universal scaling ex¬ 
ponents and scaling functions. We extracted these quan¬ 
tities for each of the momentum sectors; these are briefly 
summarized in the compilation shown in Table [ll] We 
also showed that both initial conditions containing either 
a coherent field or one with overoccupied particle modes 
arrived at the same nonthermal fixed points in each of 
the inertial ranges of momenta. 


Range 

a 

P 

7 

>^fs{pT,Pz) 

i) 

1 

2/3 

2/3 

Ai ((|p|/b)i/^ + (|p|/6)®)-i 

ii) 

-2/3 

0 

1/3 

An exp{-pll{2al))/pT 

iii) 

-1/2 

0 

1/2 

Am sech{pz/crz) 


TABLE II. Summary of the scaling properties of the three 
inertial ranges observed in the longitudinally expanding scalar 
theory, see Sec. jlllj for details. The amplitudes Ai and An 
are non-universal constants. Am is constant for Pt ^ ‘^Q 
and then rapidly decreases, and the longitudinal scale cr^ can 
depend on pr ■ 


We found that the spectrum in soft regime i) evolves 
isotropically even though the system is expanding longi¬ 
tudinally. The single particle distribution consists of a 
strong ^ |p|~^ power law enhancement. We showed that 
the distribution function in the IR follows a self-similar 
evolution. We demonstrated that our numerical results 
from the classical-statistical simulations are consistent 
with a vertex-resummed kinetic theory [69l |97] derived 
from a resummation of 2PI diagrams at NLO in an 1/N 
expansion nil Eg. Our analysis suggests that the IR dy¬ 
namics can be understood in terms of an inverse particle 
cascade from larger to lower momenta that leads to the 
formation of a Bose condensate. We found that the con¬ 
densate is formed with a power law = 0) ^ 

with the exponent a 1 in the low momentum region 
for sufficiently late times and large volumes. The con¬ 
densation time diverges with volume parametrically as 
Tc^V. Once the condensate is formed, it is subject to a 
F(p = 0) ^ power law decay. We established the 

presence of a dynamically generated mass, that decreases 
with time as m{T) ^ We found that the dynamics 

of the soft regime i) occurs for psoft < m(r); this indicates 
that the IR dynamics is that of a nonrelativistic system. 

Self-interacting scalar field theories in static geometries 
have IR dynamics reminiscent of turbulent nonrelativis¬ 
tic superfluids in 3+1-dimensions, where relativistic and 
nonrelativistic scalar field theories share a common uni¬ 
versality class [69]. The generic features of these results 
are also found for the longitudinally expanding system 
but differ from the static box results in the values of the 
dynamical exponents. 

We extended our previous analysis of the intermedi¬ 
ate momentum regime ii) [T] and verified the robustness 
of the nonthermal anisotropic fixed point by varying the 
initial conditions. For a wide range of initial conditions, 
the universal temporal scaling exponents of the longitu¬ 
dinal and transverse hard scales as well as the exponent 
characterizing the overall decrease of the occupancy are 
identical to those previously studied in depth in the longi¬ 
tudinally expanding non-Abelian gauge theory [22] . The 
functional form of the longitudinal momentum distribu¬ 
tion of the nonthermal fixed point of the scalar and gauge 
theory distributions is a Gaussian distribution and the 2- 
dimensional transverse momentum follows a thermal-like 
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behavior. This universality between the expanding 
scalar and gauge theories is truly remarkable given the 
very different scattering matrix elements in the two cases. 
Understanding this universal behavior in a microscopic 
framework remains a challenging problem. 

We discussed the emergence of a further scaling regime 
iii) in the expanding scalar theory, which appears for an 
inertial range of hard momenta Pr ^ Q- It is charac¬ 
terized by a ffattening of the single particle transverse 
momentum distribution at vanishing Pz and a hyper¬ 
bolic secant distribution along the longitudinal direction. 
The temporal evolution of the longitudinal momentum 
scale in this region of hard transverse momenta goes as 
Pz,typ ^ which is faster than the ^ scaling 

in regime ii) but distinctly slower than the ^ scal¬ 
ing expected if the system were to free stream. While 
we confirmed that this behavior corresponds to a scaling 
regime for a range of initial conditions, we are not aware 
of a known analog of this scaling in any other systems. 
Unfortunately extending the gauge theory simulations to 
the values of t/tq where this novel regime appears in 
the scalar theory simulations is not feasible with present 
computational resources. 

We also studied the implications of the fixed point 
structure for the behavior of the dynamically generated 
mass and the longitudinal and transverse pressures in 
the expanding scalar theory. We find that the domi¬ 
nant contributions to the effective mass come from the 
strongly enhanced infrared region and from the conden¬ 
sate. In contrast, in thermal equilibrium hard contri¬ 
butions at the temperature scale dominate the effective 
(Debye) mass due to the absence of a strong enhance¬ 
ment at low momenta. Similarly, soft contributions are 
important for the longitudinal pressure. We find that for 
low occupancies, where a strong infrared enhancement 
of the spectrum is not fully developed, the longitudinal 
pressure scales as as expected in kinetic theory 

from the intermediate region ii). In contrast, when the 
occupancies are large enough to generate the steep low 
momentum distribution associated with the formation of 
a condensate, we observe that the longitudinal pressure is 
distinctly different from the expectations of kinetic the¬ 
ory. This is suggestive of a strong contribution of the 
infrared sector to the longitudinal pressure. A similar 
behavior of the longitudinal pressure was seen also in the 
gauge theory case [22] ; it would be interesting to further 
explore the consequences of this finding for the validity 
of hydrodynamic descriptions. 

Our results therefore present a challenge for kinetic 
descriptions. It has been suggested that the dominance 
of small angle scattering in the expanding scalar theory 
may be because our classical-statistical simulations miss 
quantum effects that are argued to be unsuppressed for 
anisotropic systems even in the classical regime. These 
arguments further suggest that the quantum effects will 
give large contributions to the longitudinal pressure that 
are missed in the classical simulations. We showed how¬ 
ever that such quantum contributions to the longitudinal 


pressure become important precisely when the classical 
regime ends and the classical-statistical description is no 
longer valid. Hence our results are robust for / > 1 . 

A first principles description of the expanding scalar 
theory in weak coupling is desirable and should be feasi¬ 
ble. Our hope is that the universality of this dynamics 
with that of expanding non-Abelian plasmas will provide 
fundamental insight into the microscopic dynamics gov¬ 
erning the latter. 
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Appendix A: Initialization, mode functions and 
distribution function 

Within the classical-statistical description the initial 
conditions in Eqs. ^ and can be implemented by 
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sampling initial field configurations according tcp^ 


?^a(ro)+ [ J 

J 'Pt ' 


fiPT,Pz,To) + - 


¥’a(xT,?7,To) = (t)a 

X (ca.pT.i' i{PT,V,To) e*(PT^T+r^»7) + c.C.^ , (Al) 

with the abbreviations /p^ ^ = J (Ppr du/{27 t)^ and 
Pz = Here the coefficients Ca,pT,zy are Gaus¬ 

sian complex random numbers satisfying the relation 
(ca.pT.i^ 4,pT'y) = (27r)^(5a6(5(pT - Pt')<^(p - p') while 
other correlations vanish. 


Equation (Al) is a mode function decomposition of the 
scalar field (/^a(xT, tq) with mode functions 
satisfying the free equation of motion in Fourier space 


dr + ^dr +p%+'^+m^{T)] i{pT, V, t) = 0. (A2) 


Here we included for completeness a possible dynami¬ 
cally generated mass m{r), which we will neglect at ini¬ 
tial time. The mode functions are normalized as 


^(pT^y) K- C{PTy,T) = - , 


(A3) 


with adrb = a drb — b drU. With Eqs. (A2) and (AS) 


the positive frequency mode function reads 

iiPxyy) = ^ e^^/‘^Hll\pT t) , (A4) 

and its com plex conjugate is the negative frequency solu¬ 
tion of Eq. (A2). Here denotes the Hankel func¬ 


tion, which we compute numerically using the method 
described in Ref. [22]. 

Fo r v anishing transverse momentum and mass, 
Eq. (A4) is ill-defined and instead, the mode function 


can be chosen as 


«(„ = 0,.,r) = y(i) 


(A5) 


which satisfies the free equation of motion (A2) and the 
normalization condition (|A3l). 


Now let us consider later times r > tq. To compute 
the distribution function, we can use Eq. 0 with the 
statistical correlator 

Fab{PT,V,T,T') = {Pa{PT,V,T)pl{PT,V,T'))conn/V 

(A6) 

and its derivatives, where ipa{PT^i^^T) is the spatial 
Fourier transform of the classical field and V = VtL^^ 


13 yyjien considering fluctuation dominated initial conditions in 
Eq. 1^, we drop the quantum 1/2 in (Al). 


denotes the volume. Here ( )conn denotes the classical- 
statistical average of the connected correlation function. 

It is useful to get a direct relation between the statis¬ 
tical correlation function F and the distribution function 
/ that needs to be consistent with any definition of / 
including Eq. 0 . For this , we assume that a mode de¬ 
composition as in Eq. of the field is also applicable 
at later times, while interactions enter the (resummed) 
kinetic equation for the distribution function /{pr^Pzi^) 
via collision integrals. We will see in App. that this 
corresponds to a quasiparticle assumption, which is con¬ 
sistent with our considerations. For the statistical cor¬ 


relation function (A6) and its derivatives (defined above 
Eq. 0 ), the mode decomposition leads to 

F{pt,v,t) = (2/(pt,Pz,t) + 1) \i{pT,v,T)? 

F{pt,U,t) = (2f{pT,Pz,T) + 1) \Tdr^{pTy,T)f , 

(A7) 


and similar for F{pT,i^,r). With the help of Eq. (A3) 


one can easily verify that the relations (A7) are consis¬ 
tent with the definition of the distribution function for 
times r > To given by 0 . Moreover, it follows that the 
alternative definition 


fiPT,Pz,T) + ^ 


VtUN ' 


'F\\^*^PTy,T)^Pa{PTy,T) 


(A8) 


is consistent as well. The alternative definition (A8) is 


similar to what was used in the expanding gauge theory 
in Ref. [22]. We checked explicitly that both definitions 
0 and (|A8|) yield the same results at transverse mo¬ 


menta larger than the dynamically generated mass m. 
Below the mass scale, a modified dispersion should be 
taken into account. We note that our definition of the 
distribution function in Eq. 0 is independent of the spe¬ 
cific form of the mode function ^ and can thus also be 
used at low momenta. 


Appendix B: Asymptotic expressions and dispersion 
relation 


In this Appendix we consider the late-time expressions 
for the mode function '?") and statistical corre¬ 

lation function F{pT,i^,r). This provides the relation in 
Eq. (32) for finite momenta, which connects F to the dis¬ 


tribution function and to the time-dependent dispersion 
relation 


^{PT:Pz,t) = \J'm^{T)+PT + v‘^1 


(Bl) 


with (time-dependent) longitudinal momentum — vjr 
and an effective mass m(r). Similarly, the late-time ex¬ 
pressions lead to the estimate for the dispersion relation 
given by Eq. (12). 
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At late times, we make the ansatz 




1 


^2u{pt,Pz,t) 


^-iQ{pT,Pz,r) 


(B2) 


of elastic scattering reads 

d^P2 (fp3 d^Pi 




(27r)^ (27r)^ (27r)^ 


for the mode function with the time dependent function 
Vt{pT,Pz,r) = J^^dr'uj{pT,Pz,r'). To verify its validity, 
we need to show that ( |B2[ ) satisfies the free e quat ion of 
motion ( |A 2 | ) and the normalization condition ( |A3| ). The 
latte r can be easily confirmed. Since the derivatives of 
(B2) follow 


X (27r)^ + P2 - P3 - P4) 


( 5 (wi +UJ2-OJ3- OJ4) 


2i(jJ2 2C(J3 2C(J4 


(C 2 ) 


with Pi = p and with dispersion relation abbreviated 
as uji = uj{pT,i^Pz,i^^)- The effective coupling can be 
expressed as 


{d‘^ + T-'^dr)^ = -uj‘^ ^ + , (B3) Vff[/](Pl>P2,P3,P4,'r) = y [i;[/](wi +W 2 ,Pl +P 2 ,t) 


the asymptotic form in Eq. (B2) is also a solution of the 


free equation for momenta satisfying uj{pt^Pz^^) ^ 1 . 

Since in our case, the mass slowly decreases such that 

^{PTiPz^^) ^ ^ m(r)r m(ro) grows with 

time, this condition becomes eventually true for all mo¬ 
menta after a finite time. This confirms that Eq. (B2) is 


a late-time asymptotic expression for the mode function 
We no te th at alternatively, we get a simi¬ 
lar form as in Eq. (B2) by use of the tool s provided in 
Refs. j9^ llQQj when starting from Eq. (A4). 


With Eq. (B2), one finds the expressions 


\aPT,>^,r)f = 


2uj{pt,Pz,t) t 


(B4) 


and 


dT£,{PT,V,T) = -i0j{pT,Pz,T) i{pT,V,T), (B5) 


where we have dropped a term 0{^/t) in the last expres¬ 
sion. These equations, together with ( |A7|) , provide the 
estimate for the dispersion relation in ( 112 ^ and for the 


correlation function in (32) for finite momenta 


Appendix C: Collision integral from 
vertex-resummed kinetic theory 

Supplementary to the discussion of the low momen¬ 
tum region in Sec. IIIB[ we write here the expres¬ 
sion of the vertex-resummed collision integral derived in 
Refs. [69l [97] • Eollowing standard techniques [201 ESI, 
we determine its temporal scaling properties by use 
of the scaling ansatz 0 for the distribution function 

f{PT,Pz,T). 

The collision integral can be decomposed into 


C[f]ip) = J 


2-H-2 -^eff[/] ^2-0-2 

6N 


[f ], (Cl) 


with the integration measure f dO,'^^^, the effective cou¬ 
pling Agjf[/] and with the functional F^'^^[f]{{Pi}) = 
[(/i + / 2 )/ 3/4 - /i/ 2(/3 + h)]- The integration measure 


+ v[f]iLOi -W 3 ,Pi -P 3 ,r) -hu[/](wi -W 4 ,Pi - P 4 , r)] , 

(C3) 


with the abbreviation 

v[f]{0J,P,T) = 


i + n«[/](w,p,T )|2 


(C4) 


The appearance of the ‘one loop’ retarded self-energy in 
the denominator. 


n^[/](w,p,r) = 


f {pt - q.T, Pz - qz,T) 


(27r)^ 




- UJ- 


p-q 

1 


— UJ — le 


UOa — UJ- 


p-q 

1 


— UJ — le 


■ UJ- 


p-q 


T CJ T 26 CJq T UJ- 


p-q 


T CJ T 26 


(C5) 


is the result of a 1/N expansion of the 2PI effective ac¬ 
tion up to NLO and corresponds to a nonperturbative 
resummation of an infinite number of self-energy dia¬ 
grams [4^ [ 68 ] . 

The exponent p characterizing the scaling behavior of 
the collision integral C[f]{pT,Pz,r) = C[fs]{pT,Pz) 
in Eq. (22), with the abbreviations pr = r^pr and Pz = 
r^Pz , can be determined by use of the scaling ansatz ( [IT] ) 
for the distribution function. The functional 
Eq. (Cl) scales as 




impi}) = 


_ 3 q: 77 ' 2 - 0-2 


IfsKiPi}). (C 6 ) 


Eor the integration measure and the effective coupling 
we need to take into account that the dispersion rela¬ 
tion is dominated by the effective mass at low momenta, 
which leads t o Eg . ( [^ . Inserting this into the integra¬ 
tion measure (C 2 ), we have 5{uji-\-uJ2—^?,—^a) ^((Pi + 

P 2 — P 3 — P 4 j/ 2 m(r)) and 2uji 2 uj 2 2uj^ 2uj/i 16m^(r) 
to lowest nonvanishing order. According to the discus¬ 
sion below Eq. (25), we can assume an isotropic scaling 


behavior of the momenta in the IR sector with /3 = 7 . 
Then the integration measure scales as 


J dQ,‘^^‘^{{pi}) ~ T-2(2/3+7)+2/3+3a j 


(C7) 
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with the mass scaling as m{r) ^ r 


According to Eq. ( |C5| ), the retarded self-energy is 
strongly enhanced in the highly overoccupied low mo¬ 
mentum region with / ^ 1/A and therefore determines 
the scaling behavior of the effective coupling ^ 


^2 |n^| One then finds 

Collecting all parts, the collision integral scales with 
the exponent 

/i = (a — 2/3 -f cr. (C9) 
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